"""
August 21, 2015

Marine Cadastre
www.marinecadastre.gov

NOAA Office for Coastal Management
1234 South Hobson Avenue
Charleston, SC 29405
843-740-1200
ocm.mmc@noaa.gov

Software design and development by

RPS ASA, Inc.
55 Village Square Drive
Wakefield, RI 02879

IMSG
3206 Tower Oaks Blvd
Rockville, MD 20852
"""

"""
RK EDITS 5/30/2017
"""


import arcpy, os, sys, traceback, datetime, math, re

import warnings
import time
import numpy as np
from shapely.geometry import Point
from shapely.geometry.polygon import Polygon

#########################################################################################################
## Input Paramaters

##sInputFile = "" #r"E:\data\temp.gdb\day_layer2" #arcpy.GetParameterAsText(0)
##sID_Field = "MMSI" #arcpy.GetParameterAsText(1)
##sDT_Field = "BaseDateTime" #arcpy.GetParameterAsText(2)
##sOutWorkspace = "" #arcpy.GetParameterAsText(3)
##sOutName = "" # arcpy.GetParameterAsText(4)
##sLineMethod = "Maximum Time and Distance" # "Time Interval"
##sMaxTimeMinutes = 60
##sMaxDistMiles = 20
##sIntervalHours = 0
##sAddVoyageData = "false"
##voyageTable = ""
##voyageMethod = "" # "Most Frequent" "First"
##sAddVesselData = "false"
##vesselTable = ""
##
########################################
###iMaxTimeMinutes = int(sMaxTimeMinutes)
###iMaxDistMiles = int(sMaxDistMiles)
###iIntervalHours = int(sIntervalHours)
##
##prjFile = os.path.join(arcpy.GetInstallInfo()["InstallDir"],"Coordinate Systems/Geographic Coordinate Systems/World/WGS 1984.prj")
##spatialRef_WGS84 = arcpy.SpatialReference(prjFile)
##spatialRef_input = arcpy.Describe(sInputFile).spatialReference
##
###global spatialRef_CURSOR
##if spatialRef_WGS84.name <> spatialRef_input.name:
##    spatialRef_CURSOR = spatialRef_WGS84
##else:
##    spatialRef_CURSOR = None

#########################################################################################################

### Options for creating tracks
# logic for when to break line
iMaxDist = 20*0.621371       # max distance condition (miles) 

# time intervals for other conditions
iMaxTime1 = 2*60    # connect any within this # of minutes
iMaxTime2 = 1.9*60  # removing this  # connect if COG close (minutes)
iMaxTime3 = 24*60   # connect if within predicted cone (minutes) 

fSpeedMult = 0.75   # change in speed for prediction of cone
iCOGrange = 15      # COG range for prediction

iMaxCOGdiff = 15        # difference in COG for close 
iCOGdiff_predict = 25   # max COG range for matching in predicted cone 

iMaxSOGimplied = 50 # skip any points that implied SOG is very high


def RunTrackBuilder(sInputFile,sOutWorkspace,sOutName,sID_Field="MMSI",sDT_Field="BaseDateTime",sLineMethod="Maximum Time and Distance",
        iIntervalHours=0,sAddVoyageData="false",voyageTable="",voyageMethod="",sAddVesselData = "false",vesselTable = "") :
    '''
    Creating function that can be called - RK 5/30/2017
    ''' 
    arcpy.env.overwriteOutput = True
    
    #Updated to work with ArcGIS Desktop 10.1 or later
    if GetVersion() in ["10.1","10.2","10.2.1","10.2.2","10.3","10.3.1","10.4","10.4.1","10.5","10.5.1"]:
            # RK edit - added 10.5.1

        # seems to be setting a projection
        prjFile = os.path.join(arcpy.GetInstallInfo()["InstallDir"],"Coordinate Systems/Geographic Coordinate Systems/World/WGS 1984.prj")
        spatialRef_WGS84 = arcpy.SpatialReference(prjFile)
        #spatialRef_input = arcpy.Describe(sInputFile).spatialReference

        global spatialRef_CURSOR
        #if spatialRef_WGS84.name <> spatialRef_input.name:
        spatialRef_CURSOR = spatialRef_WGS84
        #else:
        #    spatialRef_CURSOR = None
        

        try:
            global bRealAISdb
            #bRealAISdb = CheckAISDatabase(sInputFile)
            bRealAISdb = True # setting to true without checking (always an AIS database)
                           
            #CreateTracks(sInputFile,sID_Field,sDT_Field,sOutWorkspace,sOutName,sAddVoyageData,sAddVesselData,sLineMethod,iMaxDistMiles,iMaxTimeMinutes,iIntervalHours)
            CreateTracks(sInputFile,sOutWorkspace,sOutName,sID_Field,sDT_Field,sLineMethod,iIntervalHours,sAddVoyageData,voyageTable,voyageMethod,sAddVesselData,vesselTable)

             
        except arcpy.ExecuteError:
            for msg in range(0, arcpy.GetMessageCount()):
                if arcpy.GetSeverity(msg) == 2:
                    arcpy.AddReturnMessage(msg)
                    print arcpy.GetMessage(msg)
         
        except:
            tb = sys.exc_info()[2]
            for i in range(0,len(traceback.format_tb(tb))):
                tbinfo = traceback.format_tb(tb)[i]
                msg = "   ERROR: \n   " + tbinfo + "   Error Info:\n   " + str(sys.exc_info()[1])

            arcpy.AddError(msg)
            print msg
    else:
        msg = "   ERROR: TrackBuilder can only be used with ArcGIS 10.1 or later"
        arcpy.AddError(msg)
        print msg
        
##########################################################################################################


def printit(inMessage):
    """Simple method to send messages to ESRI geoprocessing tool dialog, or print to screen if run from another method."""
    print inMessage
    arcpy.AddMessage(inMessage)

def CheckAISDatabase(sInputFile):
    """Performs a series of checks to see if the source points are a true AIS database with vessel and voyage data."""
    bAISdb = True
    
    desc = arcpy.Describe(sInputFile)
    input_wkspace = desc.Path
    sInName = desc.Name
     
    if os.path.splitext(input_wkspace)[1] <> ".gdb":
        bAISdb = False
    else:
        if sInName.find("Broadcast") < 0:
            bAISdb = False
        else:            
            sVesselName = re.sub("Broadcast", "Vessel", sInName)
            sVoyageName = re.sub("Broadcast", "Voyage", sInName)
    
            if arcpy.Exists(input_wkspace+"\\"+sVesselName) == False:
                bAISdb = False
            else:
                if arcpy.Exists(input_wkspace+"\\"+sVoyageName) == False:
                    bAISdb = False
                else:
                    if arcpy.Exists(input_wkspace+"\\BroadcastHasVessel") == False:
                        bAISdb = False
                    else:
                        if arcpy.Exists(input_wkspace+"\\BroadcastHasVoyage") == False:
                            bAISdb = False                         

    return bAISdb

def FormatCounts(iInCount):
    """Formats various counts to strings with comma separators, for use in dialog messages."""
    sOutCount = re.sub("(\d)(?=(\d{3})+(?!\d))", r"\1,", "%d" % iInCount)
    return sOutCount

def GetVersion():
    """Gets the version of ArcGIS being used.  Faster methods to read input datasets are used when 10.1 is the version in use."""
    dicInstallInfo = arcpy.GetInstallInfo()
    return dicInstallInfo["Version"]

def CreateDataDict_M(sInputFile,sID_Field,sDT_Field):
    """Process to read through the input point feature class.  Point information is stored in a Python dictionary in memory.
    Points are read using the Data Access Module search cursor if running ArcGIS 10.1, otherwise using a regular serach cursor if 10.0."""
    printit("\n  Reading input point data...")
    iTotPoints = int(arcpy.GetCount_management(sInputFile).getOutput(0))
    arcpy.SetProgressor("step","Reading input point data...",0,iTotPoints,1)

    dataDict = {}
    
    if bRealAISdb == True:
        #fields = ["SHAPE@XY",sID_Field,sDT_Field,"VoyageID"]
        fields = ["SHAPE@XY",sID_Field,sDT_Field,"SOG","COG","Heading","VoyageID"]

        ##FILTERING MMSI CODES THAT ARE ZERO
        where = sID_Field+" > 0"
    else:
        fields = ["SHAPE@XY",sID_Field,sDT_Field]
        where = None
     
    with arcpy.da.SearchCursor(sInputFile,fields,where,spatialRef_CURSOR,False,sql_clause=(None,'ORDER BY '+sID_Field+', '+sDT_Field+' ASC')) as rows:
        for row in rows:
            pnt = row[0]

            sID = row[1]
            dDateTime = row[2]

            if bRealAISdb == True:
                # RK Modified (correcting indices for SOG)
                iSOG = row[3]
                iCOG = row[4]
                iHeading = row[5]
                iVoyageID = row[6]
            else:
                iSOG = None
                iCOG = None
                iHeading = None
                iVoyageID = None                    
            
            if sID in dataDict:
                dataDict[sID][dDateTime] = [pnt[0],pnt[1],iSOG,iCOG,iHeading,iVoyageID]
            else:
                dataDict[sID] = {dDateTime:[pnt[0],pnt[1],iSOG,iCOG,iHeading,iVoyageID]}

            arcpy.SetProgressorPosition()
     
    printit("    Read "+FormatCounts(iTotPoints)+" points")
    arcpy.ResetProgressor()     
    return dataDict

##def CreateDataDict(sInputFile,sID_Field,sDT_Field):
##    """Process to read through the input point feature class.  Point information is stored in a Python dictionary in memory.
##    Points are read using the Data Access Module search cursor if running ArcGIS 10.1, otherwise using a regular serach cursor if 10.0."""
##    printit("\n  Reading input point data...")
##    iTotPoints = int(arcpy.GetCount_management(sInputFile).getOutput(0))
##    arcpy.SetProgressor("step","Reading input point data...",0,iTotPoints,1)
##
##    dataDict = {}
##    
##    if bRealAISdb == True:
##        fields = ["SHAPE@XY",sID_Field,sDT_Field,"VoyageID"]
##
##        ##FILTERING MMSI CODES THAT ARE ZERO
##        where = sID_Field+" > 0"
##    else:
##        fields = ["SHAPE@XY",sID_Field,sDT_Field]
##        where = None
##     
##    with arcpy.da.SearchCursor(sInputFile,fields,where,spatialRef_CURSOR,False,sql_clause=(None,'ORDER BY '+sID_Field+', '+sDT_Field+' ASC')) as rows:
##        for row in rows:
##            pnt = row[0]
##
##            sID = row[1]
##            dDateTime = row[2]
##
##            if bRealAISdb == True:
##                iVoyageID = row[3]
##            else:
##                iVoyageID = None                    
##
##            if sID in dataDict:
##                dataDict[sID][dDateTime] = [pnt[0],pnt[1],iVoyageID]
##            else:
##                dataDict[sID] = {dDateTime:[pnt[0],pnt[1],iVoyageID]}
##
##            arcpy.SetProgressorPosition()
##     
##    printit("    Read "+FormatCounts(iTotPoints)+" points")
##    arcpy.ResetProgressor()     
##    return dataDict

def CreateOuputDataset(sInputFile,sID_Field,sAddVoyageData,sAddVesselData,sOutWorkspace,sOutName):
    """Creates an empty feature class to store the track polylines.  Adds the required fields, based on the options selected by the user."""
    printit("\n  Building output track line feature class...")
    arcpy.SetProgressor("default","Building output track line feature class...")

    # RK Added to enable measure
    has_m = "ENABLED" 
     
    tmp = os.path.join('in_memory', 'template')
    #arcpy.CreateFeatureclass_management('in_memory','template',"POLYLINE","","DISABLED","DISABLED",sInputFile)
    arcpy.CreateFeatureclass_management('in_memory','template',"POLYLINE","",has_m,"DISABLED",sInputFile)

    arcpy.AddField_management(tmp,sID_Field,"LONG",)
    arcpy.AddField_management(tmp,"TrackStartTime","DATE",)
    arcpy.AddField_management(tmp,"TrackEndTime","DATE",)
     
    if bRealAISdb == True:
        if sAddVoyageData == "true":
            arcpy.AddField_management(tmp,"VoyageID","LONG",)
            arcpy.AddField_management(tmp,"Destination","TEXT")
            arcpy.AddField_management(tmp,"Cargo","LONG",)
            arcpy.AddField_management(tmp,"Draught","LONG",)
            arcpy.AddField_management(tmp,"ETA","DATE",)
            arcpy.AddField_management(tmp,"StartTime","DATE",)
            arcpy.AddField_management(tmp,"EndTime","DATE",)
        if sAddVesselData == "true":
            arcpy.AddField_management(tmp,"IMO","LONG",)
            arcpy.AddField_management(tmp,"CallSign","TEXT",)
            arcpy.AddField_management(tmp,"Name","TEXT",)
            arcpy.AddField_management(tmp,"VesselType","LONG",)
            arcpy.AddField_management(tmp,"Length","LONG",)
            arcpy.AddField_management(tmp,"Width","LONG",)
            arcpy.AddField_management(tmp,"DimensionComponents","TEXT",)

    #arcpy.CreateFeatureclass_management(sOutWorkspace,sOutName,"POLYLINE",tmp,"DISABLED","DISABLED",sInputFile)
    arcpy.CreateFeatureclass_management(sOutWorkspace,sOutName,"POLYLINE",tmp,has_m,"DISABLED",sInputFile)
     
    arcpy.Delete_management(tmp)
    del tmp
     
    arcpy.ResetProgressor()

def GetDistance(prevX,prevY,currX,currY):
    """Calculates the distance between two latitude/longitude coordinate pairs, using the Vincenty formula for ellipsoid based distances.
    Returns the distance in miles."""
    #Vincenty Formula
    #Copyright 2002-2012 Chris Veness
    #http://www.movable-type.co.uk/scripts/latlong-vincenty.html
     
    a = 6378137
    b = 6356752.314245
    f = 1/298.257223563

    L = math.radians(currX - prevX)
    U1 = math.atan((1-f)*math.tan(math.radians(prevY)))
    U2 = math.atan((1-f)*math.tan(math.radians(currY)))
    sinU1 = math.sin(U1)
    sinU2 = math.sin(U2)
    cosU1 = math.cos(U1)
    cosU2 = math.cos(U2)

    lam = L
    lamP = 9999999999
    iter_count = 0
    while abs(lam-lamP) > 1e-12 and iter_count <= 100:
        sinLam = math.sin(lam)
        cosLam = math.cos(lam)
        sinSigma = math.sqrt((cosU2*sinLam)*(cosU2*sinLam) + (cosU1*sinU2-sinU1*cosU2*cosLam)*(cosU1*sinU2-sinU1*cosU2*cosLam))

        if sinSigma == 0:
            return 0

        cosSigma = sinU1*sinU2+cosU1*cosU2*cosLam
        sigma = math.atan2(sinSigma,cosSigma)
          
        sinAlpha = cosU1*cosU2*sinLam/sinSigma
        cosSqAlpha = 1 - sinAlpha*sinAlpha
        if cosSqAlpha == 0: #catches zero division error if points on equator
            cos2SigmaM = 0
        else:
            cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha

        if cos2SigmaM == None:
            cos2SigmaM = 0

        C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha))
        lamP = lam
        lam = L+(1-C)*f*sinAlpha*(sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)))

        iter_count+=1

    uSq = cosSqAlpha*(a*a - b*b)/(b*b)
    A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)))
    B = uSq/1024*(256+uSq*(-128+uSq*(74-47*uSq)))
    deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)))
    s = b*A*(sigma-deltaSigma)

    #convert s in meters to miles
    s_miles = s*0.0006213712
    return s_miles

def GetTimeDifference(prevDT,currDT):
    """Calculates the difference between to date and time variables. The difference is returned in minutes."""
    timeDelta = currDT - prevDT
    #totMinutes = (timeDelta.days*1440) + (timeDelta.seconds/60) # original, but integer division in seconds prevents fractions of minutes
    #totMinutes = (timeDelta.days*1440) + (float( timeDelta.seconds) /60) # another option
    totMinutes = timeDelta.total_seconds() / 60 # a better option

    return totMinutes

def ProjectPoint(pnt,c,s,t) :
    # calculates prediction of next point based on cog (c) and speed(s) and time gap (t)
    # c is angle from north, s is knots; t is in hours
    # returns values in decimal degrees

    # convert c (cog) to angle from x-axis
    c90 = 90 - c  ;
    if c90<0 :
        c90 = 360 + c90  

    a = c90*math.pi / 180 ; # convert to radians

    #print a
    #print t*s*math.cos(a)
    #print t*s*math.sin(a)

    #pnt_prj  = pnt
    #print pnt_prj

    pnt_prj = [0]*2
    # dividing by 60 to convert from nm to degrees
    pnt_prj[0] = pnt[0] + ( ( t*s*math.cos(a) ) / 60 ) 
    pnt_prj[1] = pnt[1] + ( ( t*s*math.sin(a) ) / 60 ) 

    #print pnt_prj

    return pnt_prj

def degDiff(np_arr) :
    # calculating difference between two angles
    # catching and ignoring warnings
    with warnings.catch_warnings():
        warnings.simplefilter("ignore", category=RuntimeWarning)
        normDeg = np.mod( ( np_arr[1] - np_arr[0]) , 360 )
        deg_diff = min(360-normDeg, normDeg)

    return deg_diff


## RK Updated (improved logic for generating tracks, based on moore et al)
def AddLines_Segmented(dataDict,sOutputFC,sID_Field):
    """Creates track polylines and saves them in the ouput feature class.  Track lines are created using a maximum distance and
    maximum time threshold between points."""
    printit("\n  Generating track lines...")
    arcpy.SetProgressor("step","Generating track lines...",0,len(dataDict),1)
     
    cur = arcpy.da.InsertCursor(sOutputFC, ("SHAPE@",sID_Field,"TrackStartTime","TrackEndTime"))

    pnt1 = arcpy.Point()
    pnt2 = arcpy.Point()
    lineArray = arcpy.Array()
    #mArray = arcpy.Array()

    for key in sorted(dataDict.iterkeys()):
        iSegCount = 0

        #list of date time dictionary keys
        dtTimes = sorted(dataDict[key].keys())

        #start and end date variables: start date set with the start of any new line.  end date constantly updated with the second point
        #so that whenever the line ends the last datetime will already be set.
        dtSegStart = None
        dtSegEnd = None

        # get list of SOG and COG
        SOGall = [ dataDict[key][k][2] for k in dtTimes ]
        COGall = [ dataDict[key][k][3] for k in dtTimes ]

        # checking if SOG and COG are misreported (multiplied by 10)
        # SOG too high if any value greater than XX
        # COG too high if any value greater than 360 
        SOGhigh = any(i > 1020 for i in SOGall)
        COGhigh = any(i > 360 for i in COGall)

        # if either one is high then divide both COG and SOG by 10
        valMult = 1 
        if SOGhigh or COGhigh :
            valMult = 1/10 

        #print dataDict
        #print dataDict[key]
        #asdasd

        #skip MMSI if there is only one point provided
        if len(dtTimes) > 1:
            #for i in range(0,len(dtTimes)-1):
            i1 = 0# index of first point
            for i2 in range(1,len(dtTimes)):

                # getting SOG and COG as numpy arrays
                SOGvec = np.array([dataDict[key][dtTimes[i1]][2]*valMult , dataDict[key][dtTimes[i2]][2]*valMult],dtype=np.dtype('f8'))
                COGvec = np.array([dataDict[key][dtTimes[i1]][3]*valMult , dataDict[key][dtTimes[i2]][3]*valMult],dtype=np.dtype('f8'))
                HEADvec = np.array([ dataDict[key][dtTimes[i1]][4] , dataDict[key][dtTimes[i2]][4] ],dtype=np.dtype('f8'))

                # replacing with NaN if bad
                SOGvec[SOGvec == 102] = np.nan
                COGvec[COGvec == 360] = np.nan
                HEADvec[COGvec == 511] = np.nan

                pnt1.X = dataDict[key][dtTimes[i1]][0]
                pnt1.Y = dataDict[key][dtTimes[i1]][1]
                pnt1.M = time.mktime(dtTimes[i1].timetuple())
                # add M (reported SOG) if not NaN
                #if not np.isnan(SOGvec[0]) :
                #    pnt1.M = SOGvec[0]
                ## add M (time as a timestamp)
                
                pnt2.X = dataDict[key][dtTimes[i2]][0]
                pnt2.Y = dataDict[key][dtTimes[i2]][1]
                pnt2.M = time.mktime(dtTimes[i2].timetuple())
                # add M (reported SOG) if not NaN
                #if not np.isnan(SOGvec[1]) :
                #    pnt2.M = SOGvec[1]
                ## add M (time as a timestamp)

                # converting back to time
                #print datetime.datetime.fromtimestamp(pnt2.M)
                
                distance = GetDistance(pnt1.X,pnt1.Y,pnt2.X,pnt2.Y) #miles
                timeDiff = GetTimeDifference(dtTimes[i1],dtTimes[i2])
                
                SOG_implied = ( distance / (timeDiff/60) ) * 0.868976  # implied speed over ground (mph to knots)
                # getting average SOG
                # ignoring warning if both nans
                with warnings.catch_warnings():
                    warnings.simplefilter("ignore", category=RuntimeWarning)
                    SOG_reported_avg = np.nanmean(SOGvec)   # avg reported SOG
                
                #print SOG_reported_avg 

                angle_implied = math.atan2( pnt2.Y - pnt1.Y , pnt2.X - pnt1.X )  * 180 / math.pi      
                COG_implied = angle_implied + 90   ;
                if COG_implied>360 :
                    COG_implied = COG_implied - 360   

                #print COG_implied
                #print dataDict[key][dtTimes[i1]][3]
            
                # absolute value of COG difference and Heading difference
                cogDiff = degDiff(COGvec)
                headDiff = degDiff(HEADvec)
                                                    
                #SOGimplied_buf = 0.75

                # break line unless one of the following is true
                continueFlag = False

                if distance<=iMaxDist :
                    continueFlag = True
                
                elif timeDiff <= iMaxTime1 :
                    continueFlag = True

                elif timeDiff <= iMaxTime2 :
                    # if either cogDiff or headDiff is within bound
                    if cogDiff<=iMaxCOGdiff or headDiff<=iMaxCOGdiff :
                        continueFlag = True
                        
                elif timeDiff <= iMaxTime3 :
                    # calculate projected points (if a valid COG or heading)

                    # get course or heading
                    c = COGvec[0]
                    if np.isnan(c) :
                        c = HEADvec[0]

                    # only check if valid c
                    if not np.isnan(c) :
                        pnt_start = [ pnt1.X , pnt1.Y ]    
                        s = SOG_implied
                        t = timeDiff/60
                        
                        c_minus = c - iCOGrange
                        if c_minus<0 :
                            c_minus = 360 + c_minus

                        c_plus = c + iCOGrange
                        if c_plus>360:
                            c_plus = c_plus - 360

                        # Making cone
                        # speed*(1+y)
                        pnt_pred_1 = ProjectPoint(pnt_start,c_plus,s*(1+fSpeedMult),t) # cog + x 
                        pnt_pred_2 = ProjectPoint(pnt_start,c,s*(1+fSpeedMult),t) 
                        pnt_pred_3 = ProjectPoint(pnt_start,c_minus,s*(1+fSpeedMult),t) # cog - x

                        point2 = Point(pnt2.X,pnt2.Y)
                        poly_range = Polygon( [ pnt_start, pnt_pred_1, pnt_pred_2, pnt_pred_3 ] )

                        # continue if within range of predicted point
                        if poly_range.contains(point2) :
                            if cogDiff<=iCOGdiff_predict or headDiff<=iCOGdiff_predict : 
                                continueFlag = True
                            
                # skipping major anomolies (very fast implied SOG)
                skipFlag = False
                if SOG_implied > iMaxSOGimplied : 
                    skipFlag=True

##                if SOG_reported_avg >=35 : # if SOG is very high (assume SOG reported is bad)
##                    # make sure implied SOG is reasonable    
##                    if SOG_implied > 30 :
##                        skipFlag=True
##
##                else:
##                    # make sure SOG implied is close to SOG    
##                    if SOG_implied > (1+SOGimplied_buf)*SOG_reported_avg or SOG_implied < (1-SOGimplied_buf)*SOG_reported_avg :
##                        skipFlag=True

##                if continueFlag and skipFlag :
##                    print key
##                    print i1
##                    print i2
##                    print dtTimes[i1]
##                    print dtTimes[i2]
##                    print SOG_implied
##                    print SOG_reported_avg
##                    print dataDict[key][dtTimes[i2]][2]
##                    print dataDict[key][dtTimes[i1]][2]

                # update line   
                if continueFlag :
                    if not skipFlag :
                        if iSegCount == 0:
                            dtSegStart = dtTimes[i1]
                            lineArray.add(pnt1)
                            lineArray.add(pnt2)

                        else:
                            lineArray.add(pnt2)

                        dtSegEnd = dtTimes[i2]
                        iSegCount+=1

                #Break the line as the gap exceeds the thresholds
                else:
                    #print "broken"
                    if lineArray.count > 1:
                        polyline = arcpy.Polyline(lineArray,spatialRef_CURSOR,False,True)
                        cur.insertRow((polyline,key,dtSegStart,dtSegEnd))
                        del polyline

                    lineArray.removeAll()
                    dtSegStart = None
                    dtSegEnd = None

                    ##reset seg count,since line ended early due to thresholds
                    iSegCount = 0

                # iterate i1 to i2 (unless point is skipped)
                if not skipFlag :
                    i1=i2 
                
            #end line since end of this MMSI set
            if lineArray.count > 1:
                polyline = arcpy.Polyline(lineArray,spatialRef_CURSOR,False,True)
                cur.insertRow((polyline,key,dtSegStart,dtSegEnd))
                del polyline

            lineArray.removeAll()
            dtSegStart = None
            dtSegEnd = None
            #######
            
        arcpy.SetProgressorPosition()

    iTotTracks = int(arcpy.GetCount_management(sOutputFC).getOutput(0))

    printit("    Unique "+sID_Field+"'s= "+FormatCounts(len(dataDict)))
    printit("    Track lines created = "+FormatCounts(iTotTracks))
    arcpy.ResetProgressor()
    del cur


### OLD
##def AddLines_Segmented(dataDict,sOutputFC,sID_Field,iMaxDistMiles,iMaxTimeMinutes):
##    """Creates track polylines and saves them in the ouput feature class.  Track lines are created using a maximum distance and
##    maximum time threshold between points."""
##    printit("\n  Generating track lines...")
##    arcpy.SetProgressor("step","Generating track lines...",0,len(dataDict),1)
##     
##    cur = arcpy.da.InsertCursor(sOutputFC, ("SHAPE@",sID_Field,"TrackStartTime","TrackEndTime"))
##
##    pnt1 = arcpy.Point()
##    pnt2 = arcpy.Point()
##    lineArray = arcpy.Array()
##    #mArray = arcpy.Array()
##
##    for key in sorted(dataDict.iterkeys()):
##        iSegCount = 0
##
##        #list of date time dictionary keys
##        dtTimes = sorted(dataDict[key].keys())
##
##        #start and end date variables: start date set with the start of any new line.  end date constantly updated with the second point
##        #so that whenever the line ends the last datetime will already be set.
##        dtSegStart = None
##        dtSegEnd = None
##
##        #skip MMSI if there is only one point provided
##        if len(dtTimes) > 1:
##            for i in range(0,len(dtTimes)-1):
##                pnt1.X = dataDict[key][dtTimes[i]][0]
##                pnt1.Y = dataDict[key][dtTimes[i]][1]
##                pnt1.M = dataDict[key][dtTimes[i]][2] # RK - add M
##
##                pnt2.X = dataDict[key][dtTimes[i+1]][0]
##                pnt2.Y = dataDict[key][dtTimes[i+1]][1]
##                pnt2.M = dataDict[key][dtTimes[i+1]][2] # RK - add M
##
##                distance = GetDistance(pnt1.X,pnt1.Y,pnt2.X,pnt2.Y)
##                timeDiff = GetTimeDifference(dtTimes[i],dtTimes[i+1])
##
##                if distance < iMaxDistMiles and timeDiff < iMaxTimeMinutes:
##                    if iSegCount == 0:
##                        dtSegStart = dtTimes[i]
##                        lineArray.add(pnt1)
##                        lineArray.add(pnt2)
##
##                    else:
##                        lineArray.add(pnt2)
##
##                    dtSegEnd = dtTimes[i+1]
##                    iSegCount+=1
##
##                #Break the line as the gap exceeds the thresholds
##                else:
##                    #print "broken"
##                    if lineArray.count > 1:
##                        polyline = arcpy.Polyline(lineArray,spatialRef_CURSOR,False,True)
##                        cur.insertRow((polyline,key,dtSegStart,dtSegEnd))
##                        del polyline
##
##                    lineArray.removeAll()
##                    dtSegStart = None
##                    dtSegEnd = None
##
##                    ##reset seg count,since line ended early due to thresholds
##                    iSegCount = 0
##                
##                
##            #end line since end of this MMSI set
##            if lineArray.count > 1:
##                polyline = arcpy.Polyline(lineArray,spatialRef_CURSOR,False,True)
##                cur.insertRow((polyline,key,dtSegStart,dtSegEnd))
##                del polyline
##
##            lineArray.removeAll()
##            dtSegStart = None
##            dtSegEnd = None
##        arcpy.SetProgressorPosition()
##
##    iTotTracks = int(arcpy.GetCount_management(sOutputFC).getOutput(0))
##
##    printit("    Unique "+sID_Field+"'s= "+FormatCounts(len(dataDict)))
##    printit("    Track lines created = "+FormatCounts(iTotTracks))
##    arcpy.ResetProgressor()
##    del cur

def AddLines_Interval(dataDict,sOutputFC,sID_Field,iIntervalHours):
    """Creates track polylines and saves them in the ouput feature class.  Track lines are created based on the specified time interval."""
    printit("\n  Generating track lines...")
    arcpy.SetProgressor("step","Generating track lines...",0,len(dataDict),1)
     
    cur = arcpy.da.InsertCursor(sOutputFC, ("SHAPE@",sID_Field,"TrackStartTime","TrackEndTime"))

    pnt1 = arcpy.Point()
    pnt2 = arcpy.Point()
    lineArray = arcpy.Array()

    #create time interval from input parameter hours
    tdInterval = datetime.timedelta(hours=iIntervalHours)

    for key in sorted(dataDict.iterkeys()):
        #printit(key)
        iSegCount = 0

        #list of date time dictionary keys
        dtTimes = sorted(dataDict[key].keys())

        #start and end date variables: start date set with the start of any new line.  end date constatnly updated with the second point
        #so that whenever the line ends the last datetime will already be set.
        dtSegStart = None
        dtSegEnd = None

        #initialize time interval using the first point as the start, and the end based on the selected interval
        #updated as the interval threshold is reached.
        dtIntervalBegin = dtTimes[0]
        dtIntervalEnd = dtIntervalBegin + tdInterval

        #skip MMSI if there is only one point provided
        if len(dtTimes) > 1:
            for i in range(0,len(dtTimes)-1):                    
                pnt1.X = dataDict[key][dtTimes[i]][0]
                pnt1.Y = dataDict[key][dtTimes[i]][1]

                pnt2.X = dataDict[key][dtTimes[i+1]][0]
                pnt2.Y = dataDict[key][dtTimes[i+1]][1]

                if dtTimes[i] >= dtIntervalBegin and dtTimes[i+1] <= dtIntervalEnd:
                    if iSegCount == 0:
                        dtSegStart = dtTimes[i]
                        lineArray.add(pnt1)
                        lineArray.add(pnt2)
                    else:
                        lineArray.add(pnt2)

                    dtSegEnd = dtTimes[i+1]
                    iSegCount+=1

                #Break the line as the next point is outside the defined interval
                else:
                    if lineArray.count > 1:                              
                        polyline = arcpy.Polyline(lineArray,spatialRef_CURSOR)
                        cur.insertRow((polyline,key,dtSegStart,dtSegEnd))
                        del polyline
           
                    lineArray.removeAll()
                    dtSegStart = None
                    dtSegEnd = None

                    #update for next time interval
                    dtIntervalBegin = dtIntervalEnd
                    dtIntervalEnd = dtIntervalBegin + tdInterval

                    ##reset seg count,since line ended early due to thresholds
                    iSegCount = 0

            #end line since end of this MMSI set
            if lineArray.count > 1:
                polyline = arcpy.Polyline(lineArray,spatialRef_CURSOR)
                cur.insertRow((polyline,key,dtSegStart,dtSegEnd))
                del polyline
                    
            lineArray.removeAll()               
            dtSegStart = None
            dtSegEnd = None
            dtIntervalBegin = None
            dtIntervalEnd = None
        arcpy.SetProgressorPosition()

    iTotTracks = int(arcpy.GetCount_management(sOutputFC).getOutput(0))

    printit("    Unique "+sID_Field+"'s= "+FormatCounts(len(dataDict)))
    printit("    Track lines created = "+FormatCounts(iTotTracks))
    arcpy.ResetProgressor()
    del cur

def AddVoyageData_NEW(dataDict,sOutputFC,voyageTable,voyageMethod):
    """Adds voyage attribution from Voyage table to the assocatiated track lines.  Voyage data is associated with the track lines based on the VoyageID.
    Since one MMSI may have multiple VoyageIDs, the matching VoyageIDs can be selected by using the most frequent or the first."""
    printit("\n  Adding Voyage data to output...")
    
    arcpy.SetProgressor("step","Adding Voyage data to output...",0,int(arcpy.GetCount_management(sOutputFC).getOutput(0)) + int(arcpy.GetCount_management(voyageTable).getOutput(0)),1)
    
    iTotTracks = 0
    iVoyageDataAdded = 0
    
    #read ALL the voyage data into a python dictionary
    dicVoyageAttrs = {}
    with arcpy.da.SearchCursor(voyageTable,("VoyageID","Destination","Cargo","Draught","ETA","StartTime","EndTime")) as vRows:
        for vRow in vRows:
            dicVoyageAttrs[vRow[0]] = [vRow[1],vRow[2],vRow[3],vRow[4],vRow[5],vRow[6]]
            arcpy.SetProgressorPosition()
    
    #update the tracklines based on voyage data         
    with arcpy.da.UpdateCursor(sOutputFC, ("MMSI","TrackStartTime","TrackEndTime","VoyageID","Destination","Cargo","Draught","ETA","StartTime","EndTime")) as rows:
        for row in rows:
            iTotTracks+=1
               
            lMMSI = row[0]
            dtStart = row[1]
            dtEnd = row[2]

            if voyageMethod == "Most Frequent":
                #Go through data dict and create dict of voyage IDs
                dicVoyage = {}
                for dt in sorted(dataDict[lMMSI].iterkeys()):
                    if dt >= dtStart and dt <= dtEnd:
                        iVoyageID = dataDict[lMMSI][dt][2]

                        if iVoyageID in dicVoyage:
                            dicVoyage[iVoyageID]+=1
                        else:
                            dicVoyage[iVoyageID] = 1
                              
                #get most frequent Voyage of Trackline
                max_count = 0
                for v in dicVoyage.iterkeys():
                    if dicVoyage[v] > max_count:
                        max_count = dicVoyage[v]
                        iVoyageID = v

            else:#voyage method = "First"
                for dt in sorted(dataDict[lMMSI].iterkeys()):
                    if dt >= dtStart and dt <= dtEnd:
                        iVoyageID = dataDict[lMMSI][dt][2]
                        break                    

            #update row with voyageID
            row[3] = iVoyageID          

            if dicVoyageAttrs.has_key(iVoyageID):
                row[4] = dicVoyageAttrs[iVoyageID][0]
                row[5] = dicVoyageAttrs[iVoyageID][1]
                row[6] = dicVoyageAttrs[iVoyageID][2]
                row[7] = dicVoyageAttrs[iVoyageID][3]
                row[8] = dicVoyageAttrs[iVoyageID][4]
                row[9] = dicVoyageAttrs[iVoyageID][5]
    
                iVoyageDataAdded+=1

            rows.updateRow(row)

            arcpy.SetProgressorPosition()

    printit("\r    Added to "+FormatCounts(iVoyageDataAdded)+" of "+FormatCounts(iTotTracks)+" track lines")
    arcpy.ResetProgressor() 

def AddVesselData_NEW(sOutputFC,vesselTable):
    """Adds vessel attribution from Vessel table to the assocatiated track lines.  Vessel data is associated with the track lines based on the MMSI."""
    printit("\n  Adding Vessel data to output...")
    
    arcpy.SetProgressor("step","Adding Vessel data to output...",0,int(arcpy.GetCount_management(sOutputFC).getOutput(0)) + int(arcpy.GetCount_management(vesselTable).getOutput(0)),1)
    
    iTotTracks = 0
    iVesselDataAdded = 0
    
    #read ALL the vessel data into a python dictionary
    dictVesselData = {}
    with arcpy.da.SearchCursor(vesselTable,("MMSI","IMO","CallSign","Name","VesselType","Length","Width","DimensionComponents")) as vRows:
        for vRow in vRows:
            dictVesselData[vRow[0]] = [vRow[1],vRow[2],vRow[3],vRow[4],vRow[5],vRow[6],vRow[7]]
            arcpy.SetProgressorPosition()
    
    #update the tracklines based on vessel data  
    with arcpy.da.UpdateCursor(sOutputFC, ("MMSI","IMO","CallSign","Name","VesselType","Length","Width","DimensionComponents")) as rows:
        for row in rows:

            iTotTracks+=1
               
            lMMSI = row[0]
  
            if dictVesselData.has_key(lMMSI):
                row[1] = dictVesselData[lMMSI][0]
                row[2] = dictVesselData[lMMSI][1]
                row[3] = dictVesselData[lMMSI][2]
                row[4] = dictVesselData[lMMSI][3]
                row[5] = dictVesselData[lMMSI][4]
                row[6] = dictVesselData[lMMSI][5]
                row[7] = dictVesselData[lMMSI][6]

                iVesselDataAdded+=1

            rows.updateRow(row)          

            arcpy.SetProgressorPosition()

    printit("\r    Added to "+FormatCounts(iVesselDataAdded)+" of "+FormatCounts(iTotTracks)+" track lines")
    arcpy.ResetProgressor()
    
def SetupVoyageCodedDomain(sOutputFC,sInputFile,sOutWorkspace,sOutName):
    """Setup the Coded Domains for the Cargo field.  Creates a temporary table of the Cargo domain from the source database,
    and then imports it to the ouput track feature class."""
    printit("\n  Assigning Cargo coded domain values...")
    arcpy.SetProgressor("default","Assigning Cargo coded domain values...")
     
    input_wkspace = os.path.split(sInputFile)[0]

    arcpy.env.workspace = sOutWorkspace
    desc = arcpy.Describe(sOutWorkspace)
    domains = desc.domains

    #Check if domain already exists
    bCargoDomExists = False
    for domain in domains:
        if domain == "Cargo":
            bCargoDomExists = True

    if bCargoDomExists == False:
        #Copy the domain from the input AIS database to the output database
        arcpy.DomainToTable_management(input_wkspace,"Cargo","CargoDom","code","description")
        arcpy.TableToDomain_management("CargoDom","code","description",sOutWorkspace,"Cargo")

        #Delete the temporary Domain Table
        arcpy.Delete_management("CargoDom")

    #Assign domain to the field in the track feature class
    try:
        arcpy.AssignDomainToField_management(sOutName,"Cargo","Cargo")
    except:
        printit("    ERROR Assigning Domain - Workspace is read only.")
        printit("      Domain can be assigned to the 'Cargo' field manually.")
        printit("      See the 'Use Limitations' section in the help documentation for details.")

    arcpy.ResetProgressor() 


def SetupVesselCodedDomain(sOutputFC,sInputFile,sOutWorkspace,sOutName):
    """Setup the Coded Domains for the Vessel Type field.  Creates a temporary table of the VesselType domain from the source database,
    and then imports it to the ouput track feature class."""
    printit("\n  Assigning Vessel Type coded domain values...")
    arcpy.SetProgressor("default","Assigning Vessel Type coded domain values...")
     
    input_wkspace = os.path.split(sInputFile)[0]

    arcpy.env.workspace = sOutWorkspace
    desc = arcpy.Describe(sOutWorkspace)
    domains = desc.domains

    #Check if domain already exists
    bCargoDomExists = False
    for domain in domains:
        if domain == "VesselType":
            bCargoDomExists = True

    if bCargoDomExists == False:
        #Copy the domain from the input AIS database to the output database
        arcpy.DomainToTable_management(input_wkspace,"VesselType","VesselTypeDom","code","description")
        arcpy.TableToDomain_management("VesselTypeDom","code","description",sOutWorkspace,"VesselType")

        #Delete the temporary Domain Table
        arcpy.Delete_management("VesselTypeDom")

    #Assign domain to the field in the track feature class
    try:
        arcpy.AssignDomainToField_management(sOutName,"VesselType","VesselType")
    except:
        printit("    ERROR Assigning Domain - Workspace is read only.")
        printit("      Domain can be assigned to the 'VesselType' field manually.")
        printit("      See the 'Use Limitations' section in the help documentation for details.")
          
    arcpy.ResetProgressor() 
          

'''
def AddToMap(sOutputFC):
    """When running the TrackBuilder tool in ArcMap, the resulting track polylines will be loaded to the map when completed.
    Will also load the source points if they are not already loaded in the map."""
    try:
        mxd = arcpy.mapping.MapDocument("CURRENT")
        dataFrame = arcpy.mapping.ListDataFrames(mxd, "*")[0]

        #Add InputPoint if not already loaded
        desc = arcpy.Describe(sInputFile)
        sInputSource = desc.catalogPath
     
        layers = arcpy.mapping.ListLayers(dataFrame)
        bLoaded = False
        for lyr in layers:
            if lyr.dataSource == sInputSource:
                bLoaded = True
                break
            if bLoaded == False:
                lyrInput = arcpy.mapping.Layer(sInputSource)
                arcpy.mapping.AddLayer(dataFrame, lyrInput)
          
        #Add trackline
        lyrOutput = arcpy.mapping.Layer(sOutputFC)
        arcpy.mapping.AddLayer(dataFrame, lyrOutput)
    except:
        #If not run from ArcMap this will not work, and is OK
        pass
'''


#def CreateTracks(sInputFile,sID_Field,sDT_Field,sOutWorkspace,sOutName,sAddVoyageData,sAddVesselData,sLineMethod,iMaxDistMiles,iMaxTimeMinutes,iIntervalHours):
def CreateTracks(sInputFile,sOutWorkspace,sOutName,sID_Field,sDT_Field,sLineMethod,iIntervalHours,sAddVoyageData,voyageTable,voyageMethod,sAddVesselData,vesselTable):
    """Main function that calls all of the other functions.  Determines which functions are run, based on the user selected parameters."""

    #Create Dictionary of broadcast data Key = ID, Value = Dictionary of {DateTime:[x,y,voyageID]}
    #dataDict = CreateDataDict(sInputFile,sID_Field,sDT_Field)
    dataDict = CreateDataDict_M(sInputFile,sID_Field,sDT_Field)

    sOutputFC = sOutWorkspace+"\\"+sOutName
     
    #Create output feature class
    CreateOuputDataset(sInputFile,sID_Field,sAddVoyageData,sAddVesselData,sOutWorkspace,sOutName)
    #CreateOuputDataset()
          
    #Add Lines to FC
    if sLineMethod == "Maximum Time and Distance":
        #AddLines_Segmented(dataDict,sOutputFC)
        AddLines_Segmented(dataDict,sOutputFC,sID_Field) 
    elif sLineMethod == "Time Interval":
        AddLines_Interval(dataDict,sOutputFC,sID_Field,iIntervalHours)

    #Add Voyage info, if table provided:
    if sAddVoyageData == "true":
        #AddVoyageData_NEW(dataDict,sOutputFC)
        AddVoyageData_NEW(dataDict,sOutputFC,voyageTable,voyageMethod)

        #SetupVoyageCodedDomain(sOutputFC)
        #SetupVoyageCodedDomain(sOutputFC,sInputFile,sOutWorkspace,sOutName)
        

    #Add Vessel info, if table provided:
    if sAddVesselData == "true":
        #AddVesselData_NEW(sOutputFC)
        AddVesselData_NEW(sOutputFC,vesselTable)

        #SetupVesselCodedDomain(sOutputFC)
        #SetupVesselCodedDomain(sOutputFC,sInputFile,sOutWorkspace,sOutName)

    #Attempt to add layers to the map, will not work if run from ArcCatalog.
    #will not work when Run Out-of-Process
    #AddToMap(sOutputFC)
          
    printit("\n\n")


##########################################
##########################################
##if __name__ == '__main__':
##
##    arcpy.env.overwriteOutput = True
##
##    #Updated to work with ArcGIS Desktop 10.1 or later
##    if GetVersion() in ["10.1","10.2","10.2.1","10.2.2","10.3","10.3.1","10.4","10.4.1","10.5"]:
##        try:
##            global bRealAISdb
##            bRealAISdb = CheckAISDatabase()
##               
##            CreateTracks()
##             
##        except arcpy.ExecuteError:
##            for msg in range(0, arcpy.GetMessageCount()):
##                if arcpy.GetSeverity(msg) == 2:
##                    arcpy.AddReturnMessage(msg)
##                    print arcpy.GetMessage(msg)
##         
##        except:
##            tb = sys.exc_info()[2]
##            for i in range(0,len(traceback.format_tb(tb))):
##                tbinfo = traceback.format_tb(tb)[i]
##                msg = "   ERROR: \n   " + tbinfo + "   Error Info:\n   " + str(sys.exc_info()[1])
##
##            arcpy.AddError(msg)
##            print msg
##    else:
##        msg = "   ERROR: TrackBuilder can only be used with ArcGIS 10.1 or later"
##        arcpy.AddError(msg)
##        print msg
